      SUBROUTINE CONSTR(ISOLV,NCSTR,NUSO,IEND,ICY,RC,RC2,Y,AM,C1,C2,SIG,
     *Z2,Z3,C3,INCRE,FC1,FC2,C,CONLY,D)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION RC(ICY),RC2(ICY),Y(ISOLV),AM(NCSTR,IEND),C1(NCSTR,NCSTR)
     *,C2(NCSTR,NUSO),SIG(ISOLV),Z2(NCSTR),Z3(NCSTR),C3(NCSTR)
     *,C(NCSTR,IEND),CY(10),D(NCSTR,IEND)
      INTEGER TITLE(4),TITLE2(4)
C      NAMELIST/DUMP1/ RC
C      NAMELIST/INVRT/ RC
C
      LOGICAL FC1,FC2,CONLY
C
      DATA D2R/57.29577951308D0/,TITLE/'  PARTITIONED C-MATRIX  '/
C
C     STEP 1 : FORM PARTITIONED MATRIX C
      IF(FC1) GO TO 151
      DO 150 M1=1,NCSTR
      C1(M1,M1)=1.0D0
  150 CONTINUE
  151 CONTINUE
      IF(FC2) GO TO 153
      DO 152 M2=1,NCSTR
      C2(M2,M2)=-1.0D0
  152 CONTINUE
  153 CONTINUE
C
      DO 156 MR=1,NCSTR
      DO 154 MS=1,NCSTR
      C(MR,MS)=C1(MR,MS)
      D(MR,MS)=C(MR,MS)
  154 CONTINUE
      IS1=NCSTR+1
      DO 155 MS2=IS1,ISOLV
      C(MR,MS2)=C2(MR,MS2-NCSTR)
      D(MR,MS2)=C(MR,MS2)
  155 CONTINUE
  156 CONTINUE
      DO 157 MR=1,NCSTR
      C(MR,IEND)=C3(MR)
      D(MR,IEND)=C(MR,IEND)
  157 CONTINUE
C
      CALL DPMOUT(C,NCSTR,NCSTR,IEND,.FALSE.,ILAB,ILAB,24,TITLE)
C
C
C  STEP 2 : APPLY HOUSEHOLDER TRANSFORMATION TO PARTIONED MATRIX
C
      CALL THH(C,NCSTR,IEND,RC2,ICY,ISOLV,0,NCSTR,DPSOS,0)
C
COMPUTE  RC2 - INVERSE FOR THE FIRST NCSTR BY NCSTR PARAMETERS
C     STEP 3 :
      CALL UPINV(RC2,NCSTR,RC2,W)
C
C
C     MI - THE COLUMN ADDRESS IN VECTOR RC2
C
      IRC=(NCSTR*(NCSTR+1))/2
C
C     MX - THE ROW ADDRESS IN VECTOR RC AND MATRIX C2
      MA = 0
      MI = IRC
      DO 3 ICOL = 1,NUSO
      MI=MI + MA
      JJ=0
      DO 2 II=1,NCSTR
      JJ=JJ+II
      J=JJ
      C2(II,ICOL)=0.0D0
      Z2(II)=0.0D0
      DO 1 MX=II,NCSTR
      C2(II,ICOL)=C2(II,ICOL)+(-RC2(J)*RC2(MI+MX))
      Z2(II)=Z2(II)+(RC2(J)*RC2(INCRE+MX))
      J=J+MX
    1 CONTINUE
    2 CONTINUE
      MA=NCSTR+ICOL
    3 CONTINUE
C
C     FORM A NEW A-MATRIX
C     STEP 4 :
C
      DO 5 I=1,NCSTR
      DO 4 J=1,IEND
      AM(I,J)= 0.0D0
    4 CONTINUE
    5 CONTINUE
C
      MA=0
      MI=IRC
C
      DO 8 ICOL=1,NUSO
      MI=MI+MA
      JJ=0
      DO 7 II=1,NCSTR
      JJ=JJ+II
      J=JJ
      AM(II,ICOL+NCSTR)=RC(MI+II)
      Z3(II)=0.0D0
      DO 6 MX=II,NCSTR
      AM(II,ICOL+NCSTR)=AM(II,ICOL+NCSTR)+(RC(J)*C2(MX,ICOL))
      Z3(II)=Z3(II)+(RC(J)*Z2(MX))
      J=J+MX
    6 CONTINUE
    7 CONTINUE
      MA=NCSTR+ICOL
    8 CONTINUE
C
      DO 9 IZ=1,NCSTR
      AM(IZ,IEND)=RC(INCRE+IZ)-Z3(IZ)
    9 CONTINUE
C
C
C     STEP 6 : COMPUTE CONSTRAINT SOLUTION
C
      CALL THH(AM,MCSTR,IEND,RC,ICY,ISOLV,0,NCSTR,DPSOS2,0)
C      WRITE(6,DUMP1)
C
      CALL UPINV(RC,ISOLV,RC,W)
C      WRITE(6,INVRT)
C
      JJ=IRC
      ISTRT=NCSTR+1
      DO 16 IROW=ISTRT,ISOLV
      JJ=JJ+IROW
      J=JJ
      Y(IROW)=0.0D0
      SIG(IROW)=0.0D0
      DO 15 ICOL=IROW,ISOLV
      Y(IROW)=Y(IROW)+(RC(J)*RC(INCRE+ICOL))
      SIG(IROW)=SIG(IROW)+(RC(J)*RC(J))
   15 J=J+ICOL
      SIG(IROW)=DSQRT(SIG(IROW))
   16 CONTINUE
      JJ=IRC
      DO 18 IROW=1,NCSTR
      SIG(IROW)=0.0D0
      Y(IROW)=0.0D0
      DO 17 ICOL=1,NUSO
      Y(IROW)=Y(IROW)+(C2(IROW,ICOL)*Y(NCSTR+ICOL))
   17 CONTINUE
      Y(IROW)=Y(IROW)+Z2(IROW)
   18 CONTINUE
      ISTEP=0
      ISTRT=NCSTR+1
      JJ=(ISTRT*(ISTRT+1))/2
      DO 21 IROW=1,NCSTR
      SIG(IROW)=0.0D0
      ISTEP=0
      J=JJ
      DO 20 ICOL=ISTRT,ISOLV
      ISTEP=ISTEP+1
      LL=0
      DO 19 KC=1,ISTEP
      SIG(IROW)=SIG(IROW)+((C2(IROW,KC)*RC(J+LL))**2)
   19 LL=LL+1
   20 J=J+ICOL
      SIG(IROW)=DSQRT(SIG(IROW))
   21 CONTINUE
C
      IF(.NOT.CONLY) GO TO 23
      WRITE(6,25)
      DO 22 ISOUT=1,ISOLV
      WRITE(6,24) ISOUT,Y(ISOUT),ISOUT,SIG(ISOUT)
   22 CONTINUE
C
   23 CONTINUE
C
   24 FORMAT(1H ,' SOLUTION(',I4,')=',D18.9,' SIGMA(',I4,')=',D15.8)
   25 FORMAT(1H0,' PARAMETER VALUES AND SIGMAS * CONSTRAINT APPLIED')
      CALL MMAB(D,Y,NCSTR,ISOLV,1,CY,NCSTR,ISOLV,ISOLV)
      DATA TITLE2/'  CONSTRAINT  SOLUTION  '/
      CALL DPMOUT(CY,1,1,NCSTR,.FALSE.,ILAB,ILAB,24,TITLE2)
C
      RETURN
      END
